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Numerical Simulation of High-Speed 
Turbulent Reacting Flows 


F.A. Jaberi, P.J. Colucci, S. James and P. Givi 
Department of Mechanical and Aerospace Engineering 
State University of New York at Buffalo 
Buffalo, NY 14260-4400 


Abstract 

The purpose of this research is to continue our efforts in advancing the state of 
knowledge in large eddy simulation (LES) methods for computational analysis of high- 
speed reacting turbulent flows. We have just completed first year of Phase III of this 
research. This annual report provides a brief and up-to-date summary of our activities 
during the period: August 1, 1996 through July 31, 1997. 


Technical Monitor 

Dr. J. Philip Drummond (Hypersonic Propulsion Branch, NASA LaRC, Mail Stop 197, Tel: 
757-864-2298) is the Technical Monitor of this 


1 Summary of Achievements 


We are now in Year 1 of the Phase III activities on this NASA LaRC sponsored project. The 
total time allotted for this phase is three years; this phase was followed at the conclusion 
of Phase II activities (also for three years). In total we have completed 7 years of LaRC 
supported research, and two years are remaining. 

Our efforts in the past year have been particularly fruitful. We developed a new methodology 
termed the “filtered mass density function” (FMDF) and implemented it for LES of variable 
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density chemically reacting turbulent flows at low Mach numbers. The FMDF represents the 
single point joint probability density function of the subgrid scale (SGS) scalar quantities 
and is obtained by solution of its modeled transport equation. In this equation, the chemical 
reactions appear in closed form but the influences of scalar mixing and convection within the 
subgrid are modeled. The stochastic differential equations (SDEs) which yield statistically 
equivalent results to that of the FMDF transport equation axe derived and are solved via 
a Lagrangian Monte Carlo scheme. The consistency, convergence, and accuracy of FMDF 
and the Monte Carlo solution of its equivalent SDEs are assessed via comparison with data 
generated by direct numerical simulation (DNS) and with experimental data. In nonreacting 
flows, it is shown that the filtered values of temperature, density and scalars compare very 
well with those obtained by “conventional” LES procedure in which the finite difference 
solution of the transport equations governing these filtered quantities is obtained. The 
advantage of the FMDF is demonstrated in reacting flows. In the absence of closures for the 
subgrid scalar, density and temperature correlations, the results based on the conventional 
LES methods are significantly different from those based on DNS. The FMDF results show 
a much closer agreement with DNS data. The FMDF results also compare favorably with 
laboratory data of reacting turbulent shear flows, and correctly portrays several important 
features observed experimentally. 

A detailed description of the FMDF is provided in this report. 


2 Introduction 

As indicated in our last year report, (also prepared in the form of a paper by Colucci et al. 
(1997)), we have previously developed a methodology termed the “filtered density function” 
(FDF) for LES of chemically reacting flows. This methodology is based on the idea originally 
proposed by Pope (1991). The fundamental property of the FDF is to account for the effects 
of subgrid scale (SGS) scalar fluctuations in a probabilistic manner. Colucci et al. (1997) 
developed a transport equation for the FDF in which the effects of unresolved convection and 
subgrid were modeled similarly to those in conventional LES and Reynolds averaging (RA) 
procedures. This transport equation was solved numerically by a Lagrangian Monte Carlo 
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procedure and the results were compared with those obtained by direct numerical simulations 
(DNS) and by the conventional LES in different free shear flows. It was shown that in non- 
reacting flows, the first moment of the FDF as obtained by the Monte Carlo solution is the 
same as that obtained by the finite difference solution of the transport equation governing 
the mean scalar value (LES-FD). The advantage of the FDF was demonstrated in reacting 
flows in which its results were shown to deviate significantly from those based on LES-FD. 
Detailed comparison with DNS data indicated clear advantage of FDF over LES-FD. 

The encouraging results generated by FDF warrants its extensions and applications to more 
complex flows. Further assessment of its predictive capability are also in order. The objective 
in this work is to extend the methodology for treatment of reactive flows with variable density 
flows so that exothermic chemical reactions can be simulated. For that, we introduce the 
“filtered mass density function” (FMDF) which essentially is the density weighted PDF of 
SGS scalar variables. With the definition of the FMDF, the mathematical framework for 
its implementation in LES of reacting flows is established. The results obtained by FMDF 
are scrutinized by detailed comparisons with DNS and laboratory data. The FMDF deals 
only with scalar quantities; the hydrodynamic field is obtained via the conventional LES 
procedure. Probability treatment of the SGS velocity fluctuations is postponed for future 
work. Also, the formulation is based on the assumption of low Mach number. Thus while 
exothermicity and variable density effects can be studied, high speed flows cannot be treated 
by the formulation presented here. 


3 Governing Equations 

The primary transport variables in a compressible flow undergoing chemical reaction are the 
density p, the velocity vector Ui along the x; direction, the total specific enthalpy fi, the 
pressure p, and the species mass fractions Y a (a = 1,2, .. . , N 3 ). The conservation equations 
governing these variables are the continuity, momentum, enthalpy (energy) and species mass 
fraction equations, along with an equation of state relating thermodynamic variables. These 
equations are expressed as (Williams, 1985): 
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( 1 ) 
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(3) 


p = pR°T yjM Q = pUT (4) 

Of=l 

where t represents time, RP is the universal gas constant and M a is the molecular weight 
of species a. Equation (4) effectively defines the mixture gas “constant” 71. Equation (3) 
represents the mass fraction and enthalpy equations in a common form with 


4>a = Yen a = 1,2, ...,7V,, 


in which: 


N. 

4*0 = h = ^ ' h a (j) a 

a=l 



(5) 

( 6 ) 


where T denotes the temperature, To is the reference temperature and h° a and Cp a denote 
the enthalpy of formation and constant pressure specific heat of species a respectively. 


At low Mach numbers, by neglecting the viscous dissipation and thermal radiation, the 
source terms in the enthalpy equation becomes S a = Sk ^ The chemical source 

terms ( S a = S Q (<f> ), <f> = [li, V 2 , • • ■ , Yn h]) are functions of the composition variables ( <f> ). 


For a Newtonian fluid with zero bulk viscosity and Fickian diffusion, the viscous stress tensor 
r,j, mass and heat flux (J°, a = 1,2, ....,<r) are given by: 




_ ( dui 

~ ^ \d Xj 
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where fi is the dynamic viscosity and 7 = pT denotes the thermal and mass molecular 
diffusivity coefficients. In the present work, the Lewis number is assumed to be unity. 

Large eddy simulation involves the use of the spatial filtering operation (Aldama, 1990): 

/ +OO 

/(x',t)</(x',x)dx' (9) 

-OO 


where Q denotes the filter function, (f(x,t)) e represents the filtered value of the transport 
variable /(x, t), and f = f - (f) t denotes the fluctuations of / from the filtered value. For 
variable density flows it is more convenient to consider the Favre filtered quantity defined 


as: 


</(x,0)l = 


(pfh 

(p)t 


and the fluctuation about this filtered value f" = f — (f ) l . 


For a spatially and temporally invariant filter function, t/(x',x) = G(x' — x), the application 
of the filtering operation to the transport equations yields: 


dt 
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dxi 
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d(p)t{<f>a)L + d[p)i(u l )L(<l> a )L 
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dxi 


djJTU 
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(10) 

d(p)t . d(T tJ ) ( 
dxj dxi 

dTij 

dx, 

(11) 

dM a 

s' 

<* = 1)2, ,<r 

(12) 


where Tq = — (ui)i{uj)i) and M," = (p)e((ui<j> a )L — { Ui)L(<t>a)L ) denote the 

subgrid stress and the subgrid mass flux, respectively, ( pS a ) e = (p)i{S q )l (a = 1,2, ..., N„) 
are the filtered reaction rates and (pSh)e = ~ ^ e • In this work the contribution < 5/, >l is 
assumed to be negligible. 


3.1 Modeling 

In LES of non-reacting flows the closure problem is associated with Tq and M" (Erlebacher 
et al., 1992; Salvetti and Banerjee, 1995). In reacting flows, an additional model is required 
for the filtered reaction rate { S a )L . Here, modeling of (Sq)l is the subject of the probability 
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formulation as described in the next section. For the former two, we make use of currently 
available and well-established closures. The subgrid stress is modeled via two different eddy 
viscosity closures. The first is the compressible form of the Smagorinsky model (Erlebacher 
et al. y 1992): 

= — 2Cfii ({ S\})l — + gC'jl{p)<A G IIs 6 ,j (13) 

where (<S,j)i, is the resolved scale strain rate tensor, II 5 = (5;j)L(«Sij)lo A G is the charac- 
teristic size of filter and Cri, Cn are empirical constants. The drawbacks of this closure 
are well-recognized (Zang and Piomelli, 1993; Kerr et al., 1996). In an attempt to overcome 
some of these drawbacks, we also make use of a second closure which is the variable density 
form of the model that has been used in our previous work (Colucci et ai , 1997) and is given 
as: 

Tij = -2C R2 (p) ( A G E^ 2 (<«S tJ ) L - + lc 12 (p) t E6 {j (14) 

where E = \(u*)l(u*)l — ((u*)l)*' ((?/*)£,)*' |, u* = tt,— 77,- and Ui is a reference velocity in thex, 
direction. The subscript i' denotes the filter at the secondary level which has a characteristic 
size (denoted by A G ») larger than that of grid level filter. This model is essentially a modified 
version of that proposed by Bardina et al. (1983), which utilize equal sizes for the grid and 
secondary filters. We refer to this as the modified kinetic energy viscosity (MKEV) closure. 
Accordingly, the subgrid eddy viscosity, u t for these two closures are expressed as: 


= CmAyiSiMSij)!,, 

(15) 

and 


Vt = Cr2A G \J\(u*) l (u*) l - ((u*)l)/'((u*)l)**|, 

(16) 

respectively. 


A similar diffusivity model is used for the closure of the subgrid mass 

flux (Eidson, 1985): 

K = -t 

(17) 

where — (/^T*, r\ = u t /Sc t , and Sc t is the subgrid Schmidt number, assumed to be 

constant and equal to subgrid Prandtl number. 
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4 The Filtered Mass Density Function (FMDF) 


Let <p(x,t ) denote the scalar array. We define the “filtered mass density function” (FMDF), 
denoted by Fl, as: 

/ +OO 

^ / > (x , ,OC[V’,^(x , ,t)]G(x , -x)dx' (18) 


C [Vh <£(x, *)] = 6[ip - <p(x, <)] = H 6[tl> a - <p a (x, <)] ( 19 ) 

a=l 

where 8 denotes the delta function and ip denotes the composition domain of the scalar array. 
The semicolon indicates that F L (ip-,x,t) is the FMDF with regards to the compositional 
variables ip only. The term ([<£, ip(x,t)] is the “fine-grained” density (O’Brien, 1980; Pope, 
1985), and Eq. (18) implies that the FMDF is the mass weighted spatially filtered value of 
the fine-grained density. The integral property of the FMDF is such that 

/ +oo /•+ oo 

F L {ip;x, t)dip = / p(x', t)G(x! - x)dx' = (p(x, t)) t . (20) 

-OO J — oo 


The FMDF is related to the FDF (f L (ip-,x,t)) ((Colucci et ai, 1997)) through the density: 
/ , (V')/z.(V , i x ? 0 = Fi(ip] x, t). Additionally, it is useful to define the Favre-weighted FDF 
F L {rp\x,t) through (p) e F L (ip]X,t) = p(ip)f L (ip-,x,t) = F L (ip-,x,t). 

For further developments, it is useful to define the mass weighted conditional filtered mean 
of the variable Q(x,t) as: 


(Q(x,t)\ip)t 


J- ” p(x', tJQjxfi t)C [rp, <ft(x', t)] G(x' - x)dx / 

F L (ip]x,t) 


Equation (21) implies the following properties: 


( 21 ) 


(i) For Q(x,t ) = c, ( Q(x,t)\ip) t = c (22) 

(«) For Q(x,t) = Q(<p(x,t)), (Q(x,t)\ip) e = Q(ip) (23) 

f+oo 

(in) Integral property : / (Q(x,t)\ip) t F L (ip;x,t)dip 

J -OO 

= (p(x,t)) e (Q(x,t)) L (24) 
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where c is a constant, and Q(<f>(x,t)) = Q(x, t) denotes the case where the variable Q can 
be completely described by the compositional variable <f>(x, t) = [^j, <f > 2 , . . . , <f>„\. From these 
properties, it follows that the filtered value of any function of the scalar variables (such as 
p = p[</>(x, <)] and S a = S a [4>(x, <)] ) is obtained by integration over the composition space. 


To develop a transport equation for the FMDF, first the transport equation for fine-grained 
density is derived. By applying the method developed by Lundgren (1969) and others (Pope, 
1976; O’Brien, 1980) to Eq. (3) it can be shown that the fine-grained density evolves according 
to the following equation: 


dp(<f>)uj(x,t)C(4>,^) dC(<fr,V>)dJf _ d 
dt dxi dxj} Q dxi ^ difra 


[S.WCW.VO]- (25) 


The transport equation for F L {xj>\x,t) is obtained by multiplying Eq. (25) with the filter 
function G(x' - x) and integrating over x' space. The final result after some algebraic 
manipulation is 


dF L {ij>\x,t) fl[(tx,-(x,<) \if>)iF L (il>\x,t)] 

dt dxi 


d [/ i djr.. 

di pa [\ /»(<£) dxi 




dip a 


(26) 


The unclosed nature of convection and mixing is indicated by the conditional filtered value 
of the velocity and gradient of scalar flux. These unclosed terms can be further simplified. 
The convection term is decomposed into resolved and subgrid scale components as: 


{ui\i>)iF L = (ui) L F L -(- [(ui|V>)f - {ui) L \F L . (27) 

With the assumption of constant property Fickian diffusion, the conditional scalar flux gra- 
dient term may be decomposed into diffusive and dissipative parts: 

FlIp . 
(28) 

The first term on the right hand side of this equation represents the effect of molecular 
diffusion on the spatial transport of the FMDF. The second term represents the dissipative 
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d<t>* d<f > p , \ 
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nature of subgrid mixing. Upon substitution of Eqs. (27) and (28) in Eq. (26): 


9F l 

dt dxi 


jl Lmm\ _ 32 \t 

dxi \ dxi ) dip a dxl>p [ \ dxi dxi 
^[(^.1^)/ ~ (rti) L ]F L _ d[5 0 (V>)/ L ] 
dxi dxp a 


). 


IV- ) F L /p 


( 29 ) 


This is an exact transport equation for the FMDF. The last term on the right hand side 
of this equation is due to chemical reaction (and also reference pressure for a = a) and is 
in a closed form. The first term on the right hand side represents the effects of molecular 
diffusion of FDF in physical space and is closed. The second term on the left hand side 
represents convection of the FDF in physical space and is also closed provided (u,)l is 
known. The unclosed terms are associated with the second and third terms on the right 
hand side representing the influences of molecular mixing and the unresolved subgrid scale 
convection, respectively. 


The unclosed terms in Eq. (29) are modeled in a fashion consistent with conventional LES. 
The subgrid convective flux is modeled via: 

[<u.|V’)* - ( Ui) L ]F L = (30) 

The advantage of the decomposition (Eq. (27)) and the subsequent model (Eq. (30)) is 
that they yield results similar to that in conventional LES. The first two Favre moments 
corresponding to Eqs. (27) and (30) are: 


(Ui<t>a)L L ($a ) L "h L {^»)x<(^o()l]j (31) 

W/[(«i Ml - (uMMl) = -7.^^ (32) 

The term within brackets in Eq. (31) is the generalized scalar flux in the form considered 
in conventional LES (Germano, 1992; Salvetti and Banerjee, 1995). Consequently, Eq. (32) 
becomes identical to Eq. (17). 

The closure for the conditional subgrid diffusion is based on the linear mean square estimation 
(LMSE) model (O’Brien, 1980; Dopazo and O’Brien, 1976), which is also known as the IEM 
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(interaction by exchange with the mean) closure (Borghi, 1988): 


d 2 17 diadfa. A _ _ 
dxj> a dt!>i3 [ V dxi dxi ^) ( Fl/p 
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dike 


[ttm(ll> Q ~ {4 >c)l)F L ] , 


(33) 


where fi m is the “frequency of mixing within the subgrid” which is not known a priori. This 
frequency can be related to the subgrid diffusion coefficient and the filter length: fi m = 
C , n(7 + 7 t)/((^)/A^). The second moment of Eq. (33) provides an expression for the subgrid 
scalar dissipation of species a: 


£ '“ = 2 ( 7 1tef Itef), = - <*„)!) (34) 

where the subscripts in parenthesis are excluded from the summation convention. 


To establish consistency between the FMDF and conventional moment closure approaches 
we make use of the assumption: 


d_ 

dx, 



cKfVp) 

dxi 



<WM<} 

dxi 


)■ 


(35) 


Note that this assumption is not necessary to evaluate the FMDF and is only adopted to 
establish consistency. 


In the general case that the molecular diffusivity is a function of the scalar variables, the 
decomposition of the conditional scalar flux gradient takes on the form: 

Fl/p • 
(36) 

Note that the diffusive term is in a form that accounts for the correlations between the 
scalars and the molecular diffusivity. This term may be approximated by: 


_d_ 

dip o 
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1 dx i / t 


d 2 


it 


d<j>a d<j>0 


dip a dip0 [V* dxi dxi ^ 
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dwipui 

dii 
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dxi 


( p)t(T)i 
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(37) 


where (r)^, — T( {(p) l ) ■ This neglects these correlations and allows one to work directly with 
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Favre filtered quantities. This is made clear through examination of the first moment: 



d_ 

dxi 


(pW)l 


d(<t>a)_L 

dxi 


For the remainder of this work, the molecular diffusivity is assumed constant. 

With the closures given by Eqs. (30) and (33) and also the approximation made in Eq. (35), 
the modeled FMDF transport equation takes on the form: 


OFl 

dt dxi 


_d_ 

dxi 


(7 + It) 


dxi 


+ 


_d_ 

dtpc 


[n m (rp a ~ {4 >*)l)F l ) - 


d[SJF^\ 

dlp a 


(38) 


This equation may be integrated to obtain transport equations for the SGS moments. The 
equation for the first subgrid Favre moment, ( <f > a ) L , and the generalized subgrid variance, 
a l = ~ {4 >o)l are: 


0((PM4OL) + 9((p) t(uj) L{<j>a) l) 


dt 


dxi 


^(7 + 7t)~^“— + l p)t{S a )L (39) 


d((p)i°l) | ddpMuihal) 
dt dxi 


d_ ' 

dx, 


(7 + It) 



+ 2(7 + 7<) 


dj^^L d{^a)L 
dxi dxi 


2 n m (p)e<T 2 Q + 2(p)t (( <t> Q S Q ) L - (<f> Q ) L (S Q ) L ) . 


(40) 


These equations are identical to those which can be derived by filtering Eq. (3) directly, and 
adopting Eqs. (32) and (34) for the subgrid flux and dissipation. In such direct moment 
closure formulation, however, the terms involving ( S q )l remain unclosed. 


5 Lagrangian Stochastic Solution of the FMDF 

The Lagrangian Monte Carlo procedure (Pope, 1985) is employed for the solution of the 
FMDF transport equation (Eq. (38)). In this procedure, each of the Monte Carlo elements 
(particles) obeys certain equations which govern its transport. These particles undergo 
motion in physical space by convection due to the filtered mean flow velocity and diffusion 
due to molecular and subgrid diffusivities. For each particle, the compositional values of the 
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scalars are changed due to mixing and reaction. 

The spatial transport of the FMDF is represented by the general diffusion process in a 
stochastic manner using the following stochastic differential equation (SDE): (Pope, 1985; 
Risken, 1989; Gardiner, 1990) 

dXi(t) = + E(X{t),t)dWi(t ) (41) 

where X i is the Lagrangian position of a stochastic particle, Z), and E are known as the 
“drift” and “diffusion” coefficients, respectively, and W t denotes the Wiener-Levy process 
(Karlin and Taylor, 1981). The drift and diffusion coefficients are obtained by comparing 
the Fokker- Plank equation corresponding to Eq. (41) with the spatial derivative terms in the 
FMDF transport equation (Eq. (38)), 

£s\/2(7 + 7<)/M<, + (42) 

The subgrid mixing and reaction terms are implemented by altering the compositional 
makeup of the particles according to the following equation, 

^ - {Ml) + S„W + ) (43) 

where <f>+ = <f> a (Xi(t),t ) denotes the scalar value of the particle with the Lagrangian position 
vector X{. The solutions of Eqs. (41) and (43) yield the same statistics as those obtained di- 
rectly from the solution of FMDF transport equation according to the principle of equivalent 
systems (Pope, 1985; Pope, 1994). 


6 Numerical Solution Procedure 

The numerical solution of the large eddy equations is based on a hybrid procedure in which 
the hydrodynamic Favre-filtered equations (Eqs. (10) and (11)) are integrated by a finite 
difference (F.D.) method and the filtered scalar field is obtained by the Monte Carlo (M.C.) 
solution of the FMDF transport equation. The F.D. and M.C. solvers are coupled through 
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the resolved thermodynamic variables. In the presentation below, the results obtained via 
this procedure are identified by the abbreviation FMDF-1. 

For comparison, the “conventional” LES procedure is also considered. In this procedure, 
the modeled transport equations for the filtered scalar and the generalized subgrid scalar 
variance are simulated with the F.D. scheme. The hydrodynamic solver and the models for 
the subgrid stress and mass flux are identical to those in FMDF-1, but the effects of SGS 
fluctuations in the filtered reaction rate are ignored. That is, Eqs. (39)-(40) are solved via 
F.D. with the assumption {S a (^>))i = S a The results based on this procedure are 
referred to as LES-FD. 

The LES of the hydrodynamic variables, which also determines the subgrid viscosity and 
scalar diffusion coefficients, is conducted with the “compact parameter” F.D. scheme of Car- 
penter (1990). This is a variant of the MacCormack (1969) scheme in which fourth order 
compact differences are used to approximate the spatial derivatives, and a second order sym- 
metric predictor-corrector sequence is employed for time discretization. The computational 
scheme is based on a hyperbolic solver which considers a fully compressible flow. Here, the 
simulations are conducted at a low Mach number (M ss 0.3) to minimize compressibility 
effects. All the F.D. operations are conducted on fixed and equally sized grid points. Thus, 
the filtered values of the hydrodynamic variables are determined on these grid points. The 
transfer of information from these points to the location of the Monte Carlo particles is 
conducted via interpolation. Both fourth-order and second-order (bilinear) interpolations 
schemes were considered, but no significant differences in SGS statistics were observed. The 
results presented in the next section are based on simulations with fourth- and second-order 
interpolations in two-dimensional (2D) and 3D, respectively. 

The FMDF is represented by an ensemble of Monte Carlo particles, each with a set of scalars 
$jd(X( n )(<),t) and Lagrangian position vector X^. A splitting operation is employed in 
which the transport in the physical and the compositional domains are treated separately. 
The simplest means of simulating Eq. (41) is via the Euler-Maruyamma approximation 
(Kloeden and Platen, 1995): 

xntk+i) = x?(t k ) + D?(t k )At + E n (t k )(Aty^(tk) (44) 
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where D*(tk) = A(Xf n )(<*)) and E n {tk) = E(X.^(tk)). This formulation preserves the 
Markovian character of the diffusion process (Billingsly, 1979; Helfand, 1979; Gillespie, 1992) 
and facilitates affordable computations. Higher order numerical schemes for solving Eq. (41) 
are available (Kloeden and Platen, 1995), but one must be very cautious in using them 
for LES. Since the diffusion term in Eq. (41) depends on the stochastic process X(<), the 
numerical scheme must preserve the Ito-Gikhman (ltd, 1951; Gikhman and Skorokhod, 1972) 
nature of the process. Equation (44) exhibits this property. The coefficients A and E 
require the input of the filtered mean velocity and the diffusivity (molecular and subgrid 
eddy). These are provided by F.D. solution of Eqs. (10)-(11). 

The compositional values are subject to change due to subgrid mixing and chemical reaction. 
Equation (43) may be integrated numerically to simulate these effects simultaneously. Alter- 
nately, this equation is treated in a split manner. This provides an analytical expression for 
the subgrid mixing. Within a time step, At , the compositional change for the n ih particle 
due to mixing is: 

(«)*“ = + (« - («L)exp[— n m At], (45) 

Subsequently, the influence of chemical reaction is determined by evaluating the fine grain 
reaction rates and modifying the composition of the elements: 

+ At) = {4> n a ) mix + S n a At. (46) 

The implementation of Eq. (45) requires the Favre-filtered mean scalar values. These and 
other higher moments of the FMDF at a given point are estimated by consideration of 
particles within a volume centered at the point of interest. Effectively, this finite volume 
constitutes an “ensemble domain” characterized by the length scale A E (not to be confused 
with A a) in which the FMDF is discretely represented. This is necessary as, with probability 
one, no particles will coincide with the point (Pope, 1994). Here, a box of size A E is used 
to construct the statistics at the finite difference nodes. These are then interpolated to 
the particle positions. Since the mixing model only requires the input of the filtered scalar 
value, and not its derivative, this volume averaging procedure is sufficient. From a numerical 
standpoint, determination of the size of the ensemble domain is an important issue (Colucci 
et ai, 1997). Ideally, it is desired to obtain the statistics from the Monte Carlo solution 
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when the size of sample domain is infinitely small (t.e. A e — ► 0) and the number of particles 
within this domain is infinitely large. With a finite number of particles, if A# is small 
there may not be enough particles to properly construct the statistics. A larger ensemble 
domain decreases the statistical error, but increases the dispersion errors which manifests 
itself in “artificially diffused” statistical results. This compromise between the statistical 
accuracy and dispersive accuracy as pertaining to Lagrangian Monte Carlo schemes implies 
that the optimum magnitude of A^ cannot, in general, be specified a priori (Pope, 1985; 
Colucci et al., 1997). This does not diminish the capability of the procedure, but exemplifies 
the importance of. the parameters which govern the statistics. 

In an attempt to reduce computational overhead, a procedure involving the use of non- 
uniform weights is implemented. This approach allows a small number of particles to be 
imposed in regions where the behavior exhibits a low degree of variability. Conversely, in 
regions where the phenomena exhibits highly variable character, a large number of particles is 
desired to accurately represent the FMDF. This procedure is akin to grid compression in finite 
difference or finite volume approaches: whereas grid compression is a utility in increasing the 
spatial resolution in areas of interest, the application of variable weights allows for an increase 
in the resolution of the FMDF in selected regions. This is exemplified in the case of a reactive 
mixing layer. In the freestream regions where the composition is essentially described by a 
delta function, a low number of particles (in fact a single particle) is adequate to quantify 
the FMDF. In contrast, in the vicinity of the reaction zone it is desirable to impose a high 
number of particles as the composition takes on a wide range of values and the subgrid PDF 
in this region is more broad banded. In a slight generalization of that developed by Pope 
(1985), the particles evolve with a discrete Lagrangian FMDF 

W. x ;0 = Awi X) w n 8(ip - <f> n )6(x - x n ) (47) 

where Am is the mass of a particle with unit weight and w n is the weight of the n ih particle. 
Note that the semicolon is placed after the position vector x since position constitutes a 
random variable for this Lagrangian FMDF. With integration of this expression over the 
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composition domain it is possible to demonstrate: 


(P)e « 53 w * 


( 48 ) 


and the Favre filtered value of a transport quantity is constructed by the weighted 

average 


m = 


y ' («) 

Equation (48) implies that the filtered fluid density must be directly proportional to the 
sum of the weights in the ensemble domain. In the event of uniform weights, these expres- 
sions reduce to familiar forms {p)i cx Ne and (Q)l = (Pope, 1985) where Ne is 

the number of particles in the ensemble domain. Hence, with uniform weights, the particle 
number density decreases significantly in regions of high temperature. This is particularly 
important because high temperature regions are typically associated with the reaction zones 
where proper estimation of the FMDF is critical. The application of variable weights allows 
the particle number density to be increased in the high temperature reaction zones without 
having to increase the number density outside this zone. At the initial time, the number den- 
sity per grid cell of dimension A x A, (denoted NPG) is invoked. For a given value of initial 
filtered density, the particle weights can then be assigned. Similarly, at an inflow boundary, 
the weights of the particles can be assigned based upon a specified particle number density 
and mass flux. The particle weight remains fixed during its life within the computational 
domain. 


The chemical source terms in scalar equations are complex nonlinear functions of the ther- 
modynamic and compositional variables. To evaluate these source terms in the M.C. solver, 
the fine grain values of the temperature ( T n ) for all particles are calculated from the com- 
position variable <f> n = [K”, Yf , . . . , Yjf a , h n ] and the fine grain values of density (p n ) are 
determined from evaluation of the equation of state at a reference pressure po (= p n 7Z n T n ). 
The filtered pressure, ( p)t is obtained by the filtered equation of state, (p) ( = (p) t (7lT)i. 
In this equation ( p)t is obtained from the F.D. solver and the correlation ('R.T)l is obtained 
by ensemble averaging of 7 Z n T n in the M.C. solver. In this procedure, the coupling between 
the hydrodynamic and the scalar fields is taken into account and allows the investigation of 
the effects of variable density. 
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The pressure ((p)() field as determined by the above procedure exhibits some spatial os- 
cillations. It is known that in PDF methods mean quantities calculated from particles are 
accompanied by statistical noise due to the finite sample size (Pope, 1985). Since spatial 
derivatives of ( p) ( are required in the F.D. hydrodynamic solver, the oscillations can result in 
numerical difficulties. This is exacerbated by the nature of the fully compressible hydrody- 
namic code utilized in this work which allows these oscillations to propagate throughout the 
computational domain as spurious acoustic waves. Our results below show that while noise 
in the pressure field is noticeable, that of the compositional variables is not very significant. 
The amplitudes of the oscillations can be decreased by smoothing of the (' RT)l field. An 
alternate procedure is to evaluate the correlation (7 ZT)l by F.D. solution of its transport 
equation. The general form of this equation is complicated, involving two point correlations 
of scalars and temperature. At present this approach is investigated only for the case in 
which TZ is constant. With this restriction, only the solution of the Favre filtered tempera- 
ture equation is required. The reaction source term in this equation is evaluated from the 
M.C. solution. It is shown below that the statistics as determined by this procedure exhibit 
almost no spatial oscillations thus no smoothing is required. The results obtained from this 
procedure are identified by the label FMDF-2. 


7 Results 

7.1 Flows Simulated 

The FMDF is employed for the simulations of the following flow configurations: 

1. A two-dimensional (2D) temporally developing mixing layers. 

2. A 3D temporally developing mixing layer. 

3. A 2D spatially developing planar jet. 

4. A 2D spatially developing mixing layer. 

The two-dimensional simulations are conducted to allow extensive computations for assess- 
ing the consistency and accuracy of the FMDF and the convergence of the Monte Carlo 
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results. Both non-reacting and reacting flows are considered. Both the FMDF and LES-FD 
approaches are applied to the cases itemized in (l)-(3). Some of these cases are also treated 
by DNS, the results of which are used to assess the performance of the FMDF and LES-FD. 
Further appraisal is made by comparison with laboratory data for the flow under item (4). 

The temporal mixing layer consists of two co-flowing streams traveling in opposite directions 
with the same speed (Riley et al, 1986). The reactants A and B are introduced into the 
top and the bottom streams, respectively. The length in the streamwise direction is large 
enough to to allow for the rollup of two large vortices and one (subsequent) pairing of these 
vortices. In 3D simulations, the length of the domain in spanwise direction is 60% of that 
in streamwise direction (Moser and Rogers, 1991). The layer is forced via both 2D and 
3D forcing functions(Moser and Rogers, 1991; Miller et al., 1994). The initial values of the 
reactants A and B at each spanwise location in 3D simulations are identical to those in 
the 2D simulations. In the figures presented below, x, y, z correspond to the streamwise, 
cross-stream and spanwise directions, respectively. 

In the planar jet, the reactant A is injected with a high velocity from a jet of width D 
into a co-flowing stream with a lower velocity carrying reactant B (Steinberger et al., 1993). 
The size of the domain in the jet flow is 0 < x < 14£>, —3.5 D < y < 3.5 D. The ratio of 
the co-flowing stream velocity to that of the jet at the inlet is kept fixed at 0.5. A double- 
hyperbolic tangent profile is utilized to assign the velocity distribution at the inlet plane. The 
formation of the large scale coherent structures are expedited by imposing low amplitude 
perturbations at the inlet. The frequency of these perturbations correspond to the most 
unstable mode and subharmonics of this mode as determined by linear stability analysis of 
spatially evolving disturbances (Michalke, 1965; Colucci, 1993). The characteristic boundary 
condition procedure developed by Poinsot and Lele (1992) was used at the inlet. In this 
approach, the velocity and temperature variables are held fixed while the density is solved via 
its transport equation. The characteristic procedure facilitates evaluation of incoming waves 
which are necessary for the solution to the continuity equation. Zero derivative boundary 
conditions were used at the freestreams and the pressure boundary condition of Rudy and 
Strikwerda (1980) was used at the outflow. 

To demonstrate the effectiveness of the FMDF approach under more realistic conditions, 
simulations of the laboratory experiments conducted by Mungal and Dimotakis (1984) were 
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performed. In these experiments, the turbulent combustion process was studied utilizing 
a planar shear layer to sustain a diffusion flame with one freestream containing diatomic 
hydrogen and the other containing diatomic fluorine. In both streams, the reactant gas 
was diluted in nitrogen. The intention of such simulations is to verify the ability of the 
FMDF approach to handle more complex chemical kinetic mechanisms and increased values 
of the physical parameters such as the Reynolds number. While the nature of turbulence is 
inherently three-dimensional, the present simulations are restricted to two spatial dimensions. 
This assumption is not too severe since the larger, resolved scales of the flow exhibit less 
three-dimensionality than the unresolved scales (Brown and Roshko, 1974). 

In the hydrogen-fluoride mixing layer simulations, treatment of the hydrodynamic flow vari- 
ables at all of the computational boundaries is similar to that of the spatially evolving planar 
jet as described earlier. In order to mimic a naturally developing shear layer, a modified 
variant of the forcing procedure suggested by Sandham and Reynolds (1989) was utilized. 
The cross-stream velocity component at the inlet is forced at the most unstable mode as 
well as 4 harmonics (both sub- and super-) of this mode. A spatial linear stability analy- 
sis was performed to determine the most unstable mode of the hyperbolic velocity profile 
imposed at the inlet. Sandham and Reynolds (1989) used a random phase shift to ‘jitter’ 
the layer. The form of this phase shift was applied by generating a uniform random number 
—Pmax < Ay? < <pmax at each timestep and incrementing the phase by this amount. The 
effect of such a procedure is to prevent periodic behavior in the pairings of the large scale 
structures and allows the layer attain a higher degree of similarity downstream. The present 
approach considers a similar random phase shift, however, instead of utilizing a uniform dis- 
tribution for the phase shift at each time increment, a discrete approximation of the Weiner 
process was applied: 

V(t n+1 ) = ¥>(*") + y/2^At Vn , (50) 

where t/„ is a Gaussian random variable with zero mean and unit variance. The quantity a * 
is the variance of the random process per unit time and controls the ‘drift’ of the process 
over time (the variance per timestep is cr^At). 
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7.2 Reaction Mechanisms 


For the flow configurations (l)-(3), the reaction scheme is of the type A + B — * V with an 
Arrhenius expression for the reaction rate: 

S A =S B = —Aj(pA)(pB)exp(—E a /RT), (51) 

where Aj is the pre-exponential factor, E a is the activation energy, and A, B denote the 
mass fractions of species A, B , respectively. This reaction rate may be expressed in non- 
dimensional form as: 


S* A = S* B = —Da(p* A)(p* B) exp (-Ze/T*). (52) 

The asterisk denotes the non-dimensionalized variable normalized by its reference value. 
This effectively defines the Zeldovich and Damkohler numbers, Ze and Da , respectively: 


Ze = E a /RT r , 


(53) 


Da = 


Ajprfr 


(54) 


Ur/L, ' 

where p r , T r , U T and L r are reference density, temperature, velocity and length scale, re- 
spectively. Combustion exothermicity is parameterized by the non-dimensional heat release 
parameter Ce, defined as: 

-A h° P 

Ce = ^rT < 55 > 


where Ahp is the heat of formation for the product V. We note here that Ze and Ce are 
defined based on T r , the same quantity used to non-dimensionalized the temperature field 
such that T* = 1 . Other authors have chosen an alternate convention (such as the adiabatic 
flame temperature) which results in a different numerical value of Ce and Ze for the same 
flow. Both constant rate (Ze = 0) and temperature dependent (Ze ^ 0) reactions are 
considered. The species A,B,'P are assumed thermodynamically identical and the fluid is 
assumed to be calorically perfect. 


The reaction mechanism associated with the experiment is more complicated and requires 
relaxation of some of the restrictions noted in the previous paragraph. The hydrogen-fluoride 
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reaction could be represented by the reaction (Mungal and Dimotakis, 1984) 

H 2 + F 2 2 HF, A Q = 65 kcal/mol, (56) 

where A Q is the heat of reaction. The heat released in a mixture containing 1% mole fraction 
of F 2 and 1% mole fraction of H 2 diluted in nitrogen results in an adiabatic temperature rise 
of 93 1< above ambient (Mungal and Dimotakis, 1984). This reaction belongs to the more 
general family of hydrogen-halogen reactions: 

H 2 + Ha 2 — ► 2 H Ha , (57) 

in which Ha represents either F, Cl (chlorine), Br (bromine) or I (iodine). This group of 
flames has been extensively studied in the past (Chelliah, 1989) and has previously been 
used to test numerical methods for laminar flows undergoing combustion (Spalding and 
Stephenson, 1971). 

The global representation Eq. (56) is composed of a pair of second-order chain reactions 
(Mungal and Dimotakis, 1984): 

H 2 + F HF + H, AQ = 32 kcal/mol, Jfc x = 2.6 x 10 12 r°' 5 exp (58) 

BrT 

H + F 2 HF + F, AQ = 98 kcal/mol, k 2 = 3 x 10 9 T 15 exp (59) 

where the reaction rate constants and k 2 are given in units of cm 3 /(mol s), T in K , and 
the universal gas constant H° in cal/(mol K). At low concentrations of the H atom, the 
reverse of the first of these two reactions is negligibly slow (Williams, 1985). Additionally, 
the rate data suggests that the reverse of the second reaction is also negligible as compared 
to the forward reaction (Chelliah, 1989). 

The explosion limits for the hydrogen-fluorine reaction indicate that a mixture of these two 
gases at typical ambient conditions is stable (Chen et al ., 1975; Gmelin, 1980). Therefore, in 
order to initiate the reaction, a source of F atoms must be provided (Mungal and Dimotakis, 
1984). Experimentally, this was accomplished by uniformly mixing a small amount of nitric 
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oxide with the hydrogen-nitrogen mixture. The nitric oxide mixes and reacts with the fluorine 
stream to produce free fluorine atoms: 

NO + F 2 NOF + F, A Q = 18 kcal/mol, k 3 = 4.2 x 10 11 exp — (60) 

The reverse of this reaction may be neglected (Rapp and Johnston, I960). An additional 
reaction serves to limit the nitric oxide concentration (Baulch et c/., 1981; Cool et ai, 1970): 

F + NO + M NOF + M, AQ = 57 kcal/mo 1, k 4 « 3 x 10 16 cm 6 /(mol 2 s) (61) 

While it is necessary to add nitric oxide to initiate reaction, the addition of excessive amounts 
would deplete the availability of the free F atoms. Mungal and Dimotakis (1984) experi- 
mentally concluded that keeping the product of nitric oxide and diatomic fluorine molar 
concentrations to 0.03% resulted in rapid combustion. It was also noted by the authors that 
an increase of 50% in the nitric oxide concentration resulted in no appreciable changes in 
the temperature measurements. This suggests that for the conditions encountered in the 
experiment, the hydrogen-fluoride reaction can be well approximated by the limit of infinite 
rate chemistry. 

Both finite and infinite rate FMDF calculations were undertaken. For the finite rate cal- 
culations, the 8 species, 4 reaction model based on Eqs. (58)-(61) was utilized. Due to 
the very fast rate of the hydrogen-fluoride reaction, compositional change due to reaction 
was implemented in 10 chemical timesteps for every hydrodynamic timestep. The computa- 
tional cost associated with the high dimensionality as well as the large number of chemical 
timesteps motivated the investigation of infinite rate calculations. The application of infi- 
nite rate chemistry is a powerful tool in the analysis of LES flows in those cases in which 
it is justified. For a general nonlinear reaction mechanism, solution of a filtered conserved 
scalar alone is not sufficient to determine the filtered, reactive scalars since the knowledge 
of the subgrid PDF of the conserved scalar is required to determine these quantities. The 
demonstrated utility of the FMDF approach coupled with infinite rate chemistry is an added 
objective of the present work. 
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7.3 Numerical Specifications 


The results of numerical simulations are presented in three different sections. In section 7.4, 
the simulation results of the temporally developing mixing layer are used to demonstrate the 
consistency of the solution procedure for FMDF and the convergence of the Monte Carlo pro- 
cedure for the solution of FMDF. In section 7.5, the simulation results of the 3D temporally 
developing layer and those of the 2D planar jets are used to appraise the performance of the 
FMDF in predicting DNS data. In section 7.6, the simulation results of the spatially evolving 
shear layer is used to assess FMDF via comparison with experimental data. All results in 
sections 7.4 and 7.5 are presented in non-dimensional form with normalization to reference 
quantities. In all temporal flows, the reference quantities are taken to be the freestream val- 
ues, while in the spatially evolving flows, normalization is performed with regard to the high 
speed stream. The results in section 7.6 are presented in terms of dimensional quantities. In 
the presentation of the non-dimensional results the asterisks are dropped for brevity. 

The magnitude of the flow parameters considered in DNS are dictated by the resolution 
which can be afforded. The primary parameters are the flow Reynolds number (Re), the 
Damkohler number (Da), the molecular Schmidt number (5c) and the molecular Prandtl 
number ( Pr ). In all simulations Sc = Pr = 1. All finite difference simulations (in both 
DNS and LES) are conducted on equally-spaced, square (Ax = Ay — A) grids. The highest 
resolution in DNS of the 2D temporal mixing layer consists of 433 x 577 grid points which 
allows reliable calculations at Re = 11,200 (based on the velocity difference and initial 
vorticity thickness), Ce = 5, Ze = 8, and Da = 11.92. The DNS of the 3D temporal shear 
layer is conducted with a resolution of 217 x 289 x 133 with Re = 2,240, Da = 1 and 
Ce = Ze = 0. The DNS of the planar jet is performed on a 1201 x 601 grid and allows 
accurate simulations with Re = 10,000 (based on the centerline velocity at the inlet and the 
jet width), Ce = 2.5, Ze = 8 and Da = 119.2. 

The FMDF and LES-FD are conducted on lower resolution than those used in DNS. The 
LES of the temporal mixing layer is conducted on a 55 x 73 grid for 2D simulations while 
resolutions of 37 x 49 x 23 and 55 x 73 x 34 are utilized in 3D. The LES of the spatial jet 
and hydrogen-fluoride mixing layer are conducted on a 201 x 101 grid. In addition to those 
cases which are assessed via comparison with DNS, several other temporal mixing layer cases 
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are considered to demonstrate the consistency between FMDF and LES-FD approaches as 
well as convergence of the Monte-Carlo methodology. These simulations are conducted using 
the same grid resolution and Reynolds number as outlined above. A top-hat filter function 
(Aldama, 1990) of the form: 

G(x'-x) = nG(x'-Xi) 

1=1 



is used with Aq = 2A. No attempt is made to investigate the sensitivity of the results to 
the filter function (Vreman et al . , 1994) or the filter size (Erlebacher et a/., 1992). 

In both FMDF and LES-FD simulations, the subgrid stresses are modeled via the MKEV 
closure (Eq. (14)) unless specified otherwise. In the implementation of the MKEV, the 
magnitude of the reference velocity U{ is set to zero in the cross-stream and spanwise di- 
rections and to the average of the high and low speed streams in the streamwise direction. 
Additionally, the ratio of the filter size at the secondary level to that at the grid level is 
Aq'/Ag = 3. The subgrid mass flux is modeled via Eq. (17). In all cases except that 
of the hydrogen-fluoride mixing layer, Pr t = Sc t — 0.7. No attempt is made here to de- 
termine the magnitudes of the constants appearing in these models in a dynamic manner 
(Germano, 1992). In all simulations Cn = C /2 = 0.006. The magnitudes of Cm and Cm 
are given below. The subgrid mixing model requires the input of the constant Cn in the 
mixing frequency which also determines the SGS variances. A value of Cn = 4 is used in 
most simulations. Application of these values for Sc t , Pr t and Cn in the hydrogen-fluoride 
configuration resulted in a notably low temperature rise compared to the results observed 
in the experiment. To increase the subgrid diffusion and mixing, a slightly lower value of 
Set = Pr t = 0.5 and increased value of Cn = 6 were invoked. This set of constants yield 
a more favorable temperature rise. At this point, it is unclear whether these modifications 
indicate a lack of robustness in the estimation of the mixing frequency or if the low temper- 
ature rise is due to the lack of three dimensionality and small scales in the two-dimensional 
simulation. Additionally, it is noted that the mean temperature rise was found be sensitive 
to the forcing at the inlet. 
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In the FMDF of temporal mixing layer, the particles are initially distributed throughout 
the computational region. In the FMDF of the jet, the particles are supplied only in the 
region — 1.75Z? < y < 1.75 D. For the temporal shear layer, the initial value of NPG 
(defined earlier) is varied to assess its affect on statistical convergence. The size of the 
“ensemble domain” is also varied to assess its influence. The following sizes are considered: 
A e = 2A, A. The number of samples used to construct the FMDF is thus controlled by 
selection of NPG and A e- In some of the planar jet simulations and all of the hydrogen- 
fluoride mixing layer simulations, variable weights (as outlined earlier) are employed to 
decrease the overall computational expense. All other simulations utilize uniform weights. In 
the temporal mixing layer, due to flow periodicity in the streamwise and spanwise directions, 
if the particle leaves the domain at the right or the left boundary, new particles are introduced 
at the other boundary with the exact same compositional values. In the spatially evolving 
jet and mixing layer, new particles are introduced at the inlet at a rate corresponding to 
the desired (imposed) local particle number density and fluid velocity. With prescription of 
the filtered fluid density, the particle weight is adjusted to yield the proper mass flux across 
the boundary. Finally, the compositional makeup (reactant mass fractions and enthalpy) is 
imposed with a hyperbolic tangent dependence on the cross-stream coordinate. 

In application of the FMDF approach to the experimental configuration, the length scale 
used in the non-dimensionalization is taken to be 45.7 cm, corresponding to the distance 
from the virtual origin to the position downstream at which the experimental measure- 
ments were taken. The computational domain encompassed a rectangular region 1.2 by 0.6 
nondimensional units so as to decrease the effect of the outflow boundary condition at the 
(non-dimensional) position x = 1.0 coinciding with the location of the measurement appara- 
tus in the experiment. Initial discretization of the FMDF was accomplished with 5 particles 
per grid cell in the freestreams, peaking to a value of 30 at the level of the splitter plate 
(y = 0), corresponding to a initial number density from 20 to 120 sample points per ensemble 
for an ensemble domain size of A# = 2 A. The composition (reactant mass fractions and 
enthalpy) of incoming particles was set according to composition of the fluid at the point 
of entry with a hyperbolic tangent profile used to assign the inlet composition variables. 
The velocity scale was taken to be that of the fast stream (22 m/s), with a velocity ratio of 
0.4. The Reynolds number based on the downstream measurement location was taken to be 
4.0 x 10 s , consistent with the experimental value. 
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Due to the low Mach number of the experiment, very small timesteps are required to satisfy 
the CFL condition required by the compressible hydrodynamic solver. Additionally, such 
compressible solvers become sensitive to numerical instabilities at low Mach numbers. In 
order to alleviate these problems, the Mach number in the simulation was artificially in- 
creased by a factor of 5 while holding the Reynolds and Damkohler numbers fixed. This 
was accomplished by increasing the velocity reference scale by a factor of 5 and decreasing 
the reference length scale by a factor of 5. In order for the nondimensional reaction rate 
to be the same, the reaction rate constants ki through k 4 were increased by a factor of 5 2 . 
The result is an increase of the Mach number from 0.07 to 0.35 while maintaining the same 
Reynolds and Damkohler numbers. Since the experimental Mach number as well as the 
simulated Mach number are low enough such that compressibility effects on the flowfield are 
negligible, this procedure is justified. Any wave related phenomena is of course not properly 
simulated, however they have little effect on the hydrodynamics and are not of any interest 
in the present work. Dilatational effects due to heat release are not affected by this approx- 
imation and are properly accounted for. Furthermore, the filtered pressure was obtained by 
measuring the subgrid Favre correlation ( 1ZT)l from the M.C. solver and multiplying by the 
filtered density provided by the F.D. solution to the continuity equation. A minimal degree 
of smoothing of (7 ZT)l consisting of a local box filter 3x3 grid points with equal weights 
was used to minimize the magnitude of pressure oscillations throughout the flowfield. Only 
the FMDF-1 formulation was utilized in the experimental comparison. 

The simulated results are analyzed both “instantaneously” and “statistically.” In the for- 
mer, the instantaneous contours (snap-shots) and the scatter plots of the scalar values are 
considered. In the latter, the “Reynolds-averaged” statistics constructed form the instan- 
taneous data are calculated. In 2D temporal mixing layer, the flow is homogeneous in the 
x-direction; thus the statistics are constructed by the ensemble from all the grid points in 
the streamwise direction. These statistics are y dependent. In 3D shear layer the flow is 
homogeneous in x and z directions therefore statistical are evaluated by averaging in both 
these directions. In some cases, the x-averaged statistics are presented at severed different z 
locations to illustrate the three dimensional effects. In the spatially developing mixing layer 
and jet flows the averaging procedure is conducted via sampling in time. These statistics 
are y — t dependent. All Reynolds averaged results are denoted by an overbar. 
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7.4 Consistency of FMDF 


The objective in the results presented in this subsection is to demonstrate the consistency 
of the FMDF formulation. For this purpose, the LES results via FMDF and LES-FD are 
compared against each other in 2D and 3D temporal mixing layers under different conditions. 
Since the accuracy of the finite difference procedure is well-established, this comparative 
analysis provides a means of assessing the performance of the Monte Carlo solution of the 
FMDF. 

In all simulations MKEV model is used to calculate the eddy viscosity with Cm = 0.02. 
Unless otherwise mentioned, in 2D simulations NPG = 50 and in 3D simulations NPG = 20 
at locations that ( p)( = 1 . Also, in all 2D simulations A/j = A and in all 3D simulations 
A e = 2A unless specified otherwise. 

2D Simulations 

Simulations of 2D layers are conducted in which the flow is initiated with a nonuniform 
distribution of density and temperature. The initial filtered density is distributed as a 
“spike” in the middle of the layer (Fig. 1(a)). The reaction rate in this simulation is zero. 
The flow that starts with this condition evolves quite differently in comparison to one which 
starts with a uniform density distribution (Colucci, 1993). With uniform weights assigned 
to the Monte-Carlo particles, the particle number density must remain proportional to the 
fluid density. This is observed in Fig. 1(a) where it is demonstrated that the filtered density 
calculated from the Monte-Carlo particles ((p) t = ((1 /p)l )~ 1 = l/(^ £ n€A£ ^r), where 
Ne is the total number of particles involved in ensemble averaging) matches exactly with 
that of finite difference values calculated at the Eulerian grid points. Figure 1(b) shows 
that at the final time of the simulation (when the flow has experienced the pairing of two 
neighboring vortices) the Reynolds averaged filtered density calculated from finite difference 
and Monte-Carlo approaches are very close. While there are some oscillations in the particle 
number density due to finite number of particles and the stochastic nature of the Monte 
Carlo procedure, this quantity is still highly correlated with the fluid density. 

Figure 2 shows the temporal evolution of the vorticity thickness (£„) for the cases in which the 
values of the filtered density are either initially uniform or distributed as a spike. When the 
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flow starts with uniform density, the effect of thermodynamic quantities on the hydrodynam- 
ics is negligible and 8 V as obtained by FMDF-1 and FMDF-2 are identical. Comparatively, 
with an initial density spike, the growth of the layer is damped. Nevertheless, the results 
obtained by FMDF-1 are very close to those obtained by FMDF-2. The slight difference 
between the results predicted by these two methods are attributed to the differences in the 
numerical procedures. The results obtained from high resolution DNS indicate that the 
growth of layer is reasonably well predicted by the LES schemes. 

In Fig. 3, the contour plots of the resolved vorticity and temperature at the final time ( t = 44) 
as obtained by FMDF-1 and FMDF-2 are shown. Figure 3 provides a visual demonstration 
of the consistency of the FMDF as the results via the hybrid F.D. - M.C. scheme (FMDF-1) 
are in agreement with those obtained by F.D. (FMDF-2). There are some oscillations in 
the results obtained by FMDF-1. Comparatively, these oscillations are nearly eliminated in 
FMDF-2. An interesting observation in Fig. 3 is that there are significant positive vorticity 
values in the flow field. These positive values are created by the baroclinic term and are also 
observed in DNS. 

To determine the extent of noise accompanying the statistical quantities, the Reynolds av- 
eraged values of the resolved pressure, temperature, energy and conserved scalar are shown 
in Fig. 4. This figure indicates that the results based on FMDF-1 are very close to those 
obtained via FMDF-2. The most significant difference between these results is in filtered 
pressure which exhibits appreciable oscillations in FMDF-1 (Fig. 4(a)). To reduce the pres- 
sure fluctuations, a local least square filter function is applied to smooth the ( T)i values 
as obtained from the Monte Carlo particles. This operation reduces the oscillations in (p)/ 
noticeably as shown in Fig. 4(a) and does not have a significant influence on the other statis- 
tical quantities. Figures 4(b)-4(d) show that the filtered values of the resolved scale energy, 
temperature and conserved scalar as obtained by FMDF-1 are very close to those predicted 
by FMDF-2. Several other filter functions has also been employed to smooth (T)l. The 
effect of these filter functions on the flow field is demonstrated in Fig. 5, where the difference 
between calculated values of ( p) t via FMDF-2 and FMDF-1 with smoothing are shown. In 
all cases, the difference in ( p) t is very small (less than 2%). The most significant differ- 
ence between the pressure fields occurs when no smoothing operation is applied to (T)i in 
FMDF-1. The results in Fig. 5 also show that the difference between the pressure fields is 
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significantly decreased as the number of employed particles is increased. 

To demonstrate the consistency between the FMDF and moment closure (LES-FD) ap- 
proaches, comparison between the moments of a conserved scalar calculated from the Monte- 
Carlo particles and those generated by LES-FD is made. This is observed in Fig. 6 where 
contour plots of the Favre-filtered value of scalar A are presented. In the results presented in 
Fig. 6 and also in those in Figs. 7 and 8, the Favre filtered temperature in FMDF is calculated 
from Monte Carlo particles (FMDF-1) and no smoothing operation is performed. Figure 6 
clearly indicates the results obtained by FMDF are strikingly similar to those obtained by 
LES-FD. Similarly, Fig. 7 is a plot of the Reynolds averaged first and second subgrid Favre 
moments. Figure 7(a) shows the comparison of FMDF and LES-FD results for {A)l for 
different values of and NPG. The agreement in the first moment is quite good even for 
large values of A e- This is particularly important as the Favre-filtered scalar is required in 
the mixing model. The difference between the FMDF and LES-FD results is more apparent 
in Fig. 7(b) where the cross-stream variation of a\ for different A# and NPG are consid- 
ered. It is observed that FMDF overpredicts the values of the second moment. However, 
the difference between FMDF and LES-FD diminishes as A# decreases. This is consistent 
with the results obtained in constant density flows (Colucci et al., 1997). 

The consistency between FMDF and LES-FD as established above is only valid in nonre- 
acting flows. In reacting flows, the reaction term in FMDF transport equation is closed but 
the corresponding terms in the moment equations require modeling. This lack of consis- 
tency is demonstrated in Fig. 8 where simulations of a uniform density reacting temporal 
mixing layer with Da — 2 and Ce = Ze = 0 are considered. In these simulations, the LES 
resolution is lowered to 37 x 49 and Re = 2, 800. The product thickness (<5 P ) from FMDF 
and Favre-filtered DNS are nearly identical. Comparatively, in the absence of a reaction 
closure, the product formation in LES-FD is significantly higher. Also shown in Fig. 8 are 
the results calculated from the filtered density function (FDF) (Colucci et al., 1997). The 
results obtained from FMDF and FDF are quite close, indicating the consistency of FMDF 
with FDF formulation in the limit of zero density variations. 

3D Simulations 

To further demonstrate the consistency between the FMDF and LES-FD methodologies, 
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large eddy simulations of 3D temporal shear layers are conducted. In these simulations the 
initial density and temperature profiles are distributed as a “spike”, similar to that in the 
previous 2D simulations. 

The results of 3D simulations with 2D forcing are similar to those obtained via 2D simu- 
lations. To demonstrate this similarity, in Fig. 9 the cross-stream variation of the filtered 
density averaged over the x direction at several different z locations as obtained by FMDF-2 
is shown. The values are calculated from the Monte- Carlo particles (M.C.) and finite differ- 
ence procedure (F.D.). Figure 9 demonstrates that the filtered density profiles at different 
spanwise locations are similar to each other and to the results shown in Fig. 1(b) for 2D 
simulations. All other statistics were observed to exhibit similar two-dimensional behavior. 
Consistent with the 2D simulations, Fig. 9 further indicates that there is an excellent agree- 
ment between the filtered density calculated from M.C. and F.D. The agreement obtained 
by the different methods is better illustrated in Fig. 10, where scatter plots of ( T) L and ( A) L 
as calculated by M.C. and F.D. are shown. It is observed that the local values of ( T) L and 
(A)l calculated from Monte Carlo particles are highly correlated with those calculated by 
the finite difference method. The correlation coefficients among the data presented in Figs. 
10(a) and 10(b) are 0.998 and 0.993, respectively. 

The statistics calculated from simulations with 3D forcing exhibit significant variations along 
the spanwise direction. This is observed in Fig. 11, where the cross-stream variation of the 
x-averaged streamwise vorticity and pressure at several z locations obtained by FMDF-1 and 
FMDF-2 are shown. The results in Fig. 11(a) indicate that the filtered vorticity calculated 
by FMDF-1 is very close to that calculated by FMDF-2. The filtered pressure obtained from 
FMDF-1 also exhibits similar trends to that obtained from FMDF-2, however, the noise in 
the pressure field as calculated by FMDF-1 is appreciable. Despite the obvious effect of this 
noise on the pressure field, the filtered scalar and temperature values as calculated by M.C. 
are consistent with those calculated by F.D. This is illustrated Fig. 12 in which scatter plots 
of ( T)l and ( A)l are shown. The results in Fig. 12(a) indicate that the local values of ( T)l 
calculated from Monte Carlo particles are very close to those calculated by finite difference 
method. Furthermore, there is a high correlation between the local values of {A) L calculated 
from the Monte-Carlo particles and finite difference (Fig. 12(b)). The correlation coefficients 
among the data presented in Figs. 12(a) and 12(b) are 0.999 and 0.999, respectively. 
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7.5 Validations via DNS 


The objective of this section is to assess the overall performance of FMDF methodology, to 
appraise the validity of the submodels employed in the FMDF transport equation, and to 
demonstrate the capabilities of FMDF for LES of complex exothermic chemically reacting 
turbulent flows. For this objective, the FMDF results are compared against DNS of the 
same flow configuration with the same magnitudes of the physical parameters (Re, Da 
etc...). For meaningful comparison, the DNS data is filtered and down-sampled onto coarse 
grids corresponding to the resolution used in FMDF. At this point it should be emphasized 
that FMDF is not intended to be an alternative to DNS, rather the comparison is made to 
assess the performance of the FMDF approach. The value in the FMDF approach lies in 
its applicability to more “practical” systems which are incapable to be treated by DNS. To 
illustrate the capability of the FMDF, the results are also compared with LES-FD in which 
the effects of SGS fluctuations on the filtered reaction rate are ignored. The results of two- 
and three-dimensional simulations are discussed in different subsections. First the results 
for the 2D temporal shear layer and 2D planar jet are presented. This is followed by the 3D 
temporal shear layer results. Unless otherwise specified, all FMDF calculations presented in 
this section utilize the FMDF-2 formulation. 

2D Simulations 

To quantify the performance of the FMDF methodology in the case of heat release, simula- 
tions of a temporal mixing layer undergoing exothermic chemical reaction are conducted. In 
Fig. 13(a), the FMDF predictions of the product thickness are compared with DNS results. 
The cross-stream variation of the product mass fraction averaged over streamwise direction 
at t = 44 are shown in Fig. 13(b). The FMDF results are calculated with both A*; = A 
and A e = 2A. The results obtained by LES-FD are also presented to demonstrate the 
importance of the SGS scalar fluctuations. Figure 13(a) indicates that the values of the 
product thickness predicted by FMDF are reasonably close to those of DNS. The size of the 
ensemble domain used to calculate the Favre-filtered statistics from particles does not have 
significant influence on the magnitude of the product thickness. Additionally, the results in 
Fig. 13(a) illustrate that the neglect of SGS correlations results in significant overpredictions 
of the product mass fraction. This behavior is observed at all times and and is consistent 
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with that in Reynolds averaging (Libby and Williams, 1980). The results shown in Fig. 
13(b) further illustrate these trends. While the cross-stream product distribution obtained 
in DNS and FMDF with either Aj g = A or Ae = 2A are in good agreement, that predicted 
by LES-FD is significantly higher. It is expected that the deviation between DNS and LES- 
FD results would increase as the magnitude of the Damkohler number and/or Reynolds 
number increases (Colucci et ah, 1997). In the majority of “practical” reacting systems, 
the magnitudes of the Damkohler and Reynolds numbers are considerably large. Therefore 
it is expected that the effects of the SGS correlations would be highly pronounced in such 
applications. 

In combusting flows, the reaction rate is typically a strong function of temperature. Therefore 
accurate prediction of the scalar field relies on the proper assessment of the temperature. In 
Fig. 14 the cross-stream variation of the filtered temperature (averaged over the x direction) 
for DNS, FMDF and LES-FD at t = 44 are shown. This figure shows that for this flow with 
significant variations in temperature due to heat release, the averaged filtered temperature is 
reasonably well predicted by FMDF. In contrast, LES-FD yields a significant overprediction 
in the product formation, indicating the importance of the SGS scalar and temperature 
fluctuations. Furthermore, it is demonstrated in Fig. 14 that the filtered temperature as 
obtained from M.C. particles is nearly identical with that obtained by the F.D. solution. 

Accurate prediction of the filtered temperature field is also critical in evaluation of the 
velocity field. In the cases considered here the heat release associated with the reaction has 
a significant influence on the hydrodynamics. This is observed in Fig. 15 where it is shown 
that the vorticity thickness (6 V ) calculated by DNS is significantly decreased due to heat 
release. Figure 15 also shows that the vorticity thickness obtained from FMDF exhibits 
similar trends to those of DNS. Nevertheless, in both heat releasing and non-heat releasing 
cases, FMDF underpredicts the DNS values due to inaccuracy of the models for SGS stresses. 
The inaccuracy of the SGS stress models has a noticeable effect on the development of the flow 
in the transitional period of its evolution and is the primary reason for the slight differences 
between the DNS and FMDF results in Figs. 13 and 14. 

The effectiveness of the FMDF to predict the slightly more complex jet flow is summarized 
in Figs. 16-21. As described above, in the treatment of exothermic flows, 3 different method- 
ologies for dealing with the filtered pressure are employed. In each procedure, the filtered 
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density obtained from the finite difference solution of the filtered equations is used in the 
equation of state. In the first approach, the Favre-filtered temperature obtained by ensemble 
averaging over the Monte-Carlo particles is used to evaluate the filtered pressure (FMDF- 
1). Due to the statistical noise of the finite samples utilized to calculate mean quantities, 
the resulting pressure field suffers from oscillations. This is demonstrated by a snapshot of 
the flowfield in Figs. 16(a) and 17(a) in which contours of the pressure and temperature 
for exothermic reacting simulations are displayed. This behavior is amplified by the fact 
that the resulting pressure field is differentiated in the momentum equations and the finite 
difference solver is acoustic in nature. In an attempt to reduce this oscillatory behavior, a 
local box filter is applied to smooth the temperature. Figure 16(b) demonstrates that the 
pressure fluctuations in the resulting flowfield are reduced significantly. Furthermore, addi- 
tional diffusion due to the smoothing nature of the box filter is apparent, especially in the 
temperature contours (Fig. 17(b)). The third approach considers the evaluation of the Favre 
filtered temperature by considering its transport equation (FMDF-2). The subgrid corre- 
lations appearing as the source term in this equation are obtained from the FMDF solver. 
The resulting pressure field, depicted in Fig. 16(c), is nearly free of pressure oscillations. 
The temperature values evaluated via this procedure are also free of noise (Fig. 17(c)). It 
is to be noted that in the present simulations, the gas constant of the mixture is assumed 
constant and any mass fraction-temperature correlations in the filtered equation of state are 
identically zero. The numerical treatment of these correlations due to variable mixture gas 
constant is the focus of our current investigation. 

It has been shown in Fig. 1 that in the nonreacting density-stratified mixing layer, the 
Reynolds averaged values of the particle number density correlate with the filtered fluid 
density. Similarly, Fig. 18 shows that in reacting planar jet simulations the instantaneous 
particle number density and the filtered fluid density as calculated by FMDF are reasonably 
correlated. The results considered in Figs. 1 and 18 correspond to cases in which the particles 
have uniform weights. Of particular interest in Fig. 18 is that the particle number density 
is lowest in the high temperature reaction zones. To assess the effectiveness and accuracy 
of the variable weight procedure, simulations of a reacting planar jet identical to that in 
Figs. 16-18 but with variable particle weights are conducted. Particles are introduced at 
the inlet based on a specified number density which is lowest in the freestreams and peaks 
at the reaction zones. The resulting particle number density, sum of particle weights, and 
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filtered fluid density calculated by FMDF are shown in Fig. 19. It is observed in Fig. 19(a) 
that downstream there is a higher proportion of particles in the reaction zones relative to 
the case with uniform weights. The sum of the particle weights as shown in Fig. 19(b) is 
reasonably correlated with the filtered fluid density (Fig. 19(c)). A comparison between 
Figs. 18(b) and 19(c) indicate that despite the significant difference in the total number of 
particles and particle weighing procedures, the filtered density fields are nearly identical in 
the two simulations. The computational time in the case with variable weights is lower than 
that in the case with equal particle weighting by nearly a factor of two. 

The similarity of the results obtained with uniform and variable particle weighting procedures 
is further demonstrated in Fig. 20(a) where the variation of the integrated product thickness 
with the downstream coordinate is shown. Also shown in this figure is the product thickness 
as obtained from DNS and LES-FD. It is observed that LES-FD substantially overpredicts 
the amount of product in comparison to the filtered DNS data. The FMDF results are much 
closer to those of DNS. Further comparison between DNS, FMDF and LES-FD results is 
made Fig. 20(b) where the time-averaged profile of the product variation in y direction at 
x/D = 14 is shown. Similar to the conclusions drawn from Fig. 20(a), it is observed that 
the amount of product predicted via LES-FD is significantly higher than that of DNS while 
the FMDF results are close to DNS values. 

The ability of the FMDF to capture the effects of chemical reaction lies in its ability to resolve 
the scalar and temperature correlations without modeling. These correlations constitute an 
important contribution to the filtered reaction rate. To demonstrate this, in Fig. 21, contour 
plots of the SGS “unmixedness”, defined as (£(<£)), - S((<f>) L ), from both DNS and FMDF 
are shown. A comparison between the results shown in this figure clearly indicate that 
the SGS unmixedness calculated by FMDF is in very good agreement with that of DNS. 
Comparatively, this term is effectively ignored in LES-FD. The contribution of the SGS 
unmixedness to the total filtered reaction rate is expected to increase as the flow Reynolds 
and Damkohler numbers increases. Therefore, it is expected that the difference between 
DNS and LES-FD results increase as Reynolds and/or Damkohler numbers increase. 

3D Simulations 

The major conclusions drawn from the 2D results axe confirmed in 3D simulations. To 
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build upon the results obtained from the two-dimensional simulations, DNS and LES of 
three-dimensional temporal shear layer flow are conducted. The temporal evolution of the 
product thickness as predicted by FMDF is compared with DNS and LES-FD results in Fig. 
22 for both low and high LES resolutions. Figure 22(a) illustrates the importance of the 
SGS unmixedness in 3D flows, as the product thickness obtained from LES-FD is appreciably 
higher than that predicted by DNS. Consistent with the results of the 2D simulations, the 
product thickness predicted by FMDF with MKEV model are close to those of DNS. Also 
shown in Fig. 22(a) are the FMDF results obtained using the Smagorinsky model as the 
SGS stress closure. The accuracy of the Smagorinsky model for this transitional flow is 
less than that of the MKEV model as the product thickness values obtained by FMDF 
with Smagorinsky closure compare less favorably to the DNS results than those obtained by 
FMDF with the MKEV model. The results shown in Fig. 22(b) for higher LES resolutions 
resemble those shown in Fig. 22(a) for lower LES resolution. Again LES-FD significantly 
overpredicts the DNS results while the results obtained from FMDF using the MKEV model 
yields a more favorable comparison. Expectedly, the difference among DNS, LES-FD and 
FMDF results decrease as the resolution in LES increases. 

To further demonstrate the capabilities of FMDF, in Fig. 23 the cross-stream variations of 
the filtered product mass fractions for DNS, LES-FD and FMDF at t = 44 are shown. The 
product values reported in Fig. 23(a) are obtained by averaging the instantaneous results 
over x and 2 directions. To illustrate the three-dimensionality of the flow, the product mass 
fractions averaged over x direction at two different 2 locations are presented in Fig. 23(b). 
Consistent with the results shown in Fig. 22, it is observed that the product mass fractions 
are significantly overpredicted by LES-FD, again indicating the importance of the SGS scalar 
correlations. The product mass fraction predicted by FMDF as shown in Fig. 23(a) are in 
agreement with that of DNS. The FMDF reasonably predicts the product mass fractions for 
several spanwise locations demonstrating the enormous potential of the FMDF for LES of 
three-dimensional complex turbulent reacting flows. 

A common feature of the results shown in Fig. 23 is that the dispersion of the product 
values along the cross-stream direction is slightly underpredicted by FMDF. The reason for 
this behavior is that the SGS stress and scalar flux models are not entirely accurate and 
result in underprediction of the growth of the layer. This is illustrated in Fig. 24, where it 
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is shown that for different filter sizes (different LES resolution) the temporal evolution of 
the vorticity thickness (<$„) as predicted by FMDF with the MKEV model are lower than 
that that predicted by DNS. Expectedly, as the resolution of FMDF is increased, the results 
better compare with DNS. The vorticity thickness predicted by FMDF using the Smagorinsky 
model are also considered in Fig. 24. It is shown that the vorticity thickness predicted by 
FMDF with the MKEV model predicts the DNS data much better than that via FMDF 
with the Smagorinsky closure. It is well established that the Smagorinsky closure generates 
excessive damping on the resolved scales in transitional regions (Zang and Piomelli, 1993) 
which, consequently, yields wrong prediction of the growth rate of the shear layer. The 
improved prediction of the eddy viscosity also results in more accurate prediction of the 
product formation by FMDF as illustrated in Fig. 22. This is primarily due to the role of 
the eddy viscosity in the determination of the subgrid mixing frequency. 

Computational Requirements 

To evaluate the computational cost of FMDF, the computational times associated with the 
simulations of 2D reacting planar jet and 3D reacting mixing layer are shown in tables 1 and 
2, respectively. The overhead of the FMDF simulation is somewhat extensive as compared 
to LES-FD; nevertheless the running time for FMDF simulation is significantly less than 
that of DNS. While this overhead is tolerated in present simulations, there are several means 
of reducing it for future applications. For example the Lagrangian procedure would benefit 
from the utilization of parallel architecture, since a significant portion of the time is devoted 
to computations in large loops dimensioned by the total number of Monte Carlo particles. 
This has been discussed for use in PDF methods(Leonard and Dai, 1994) and its utilization 
in FMDF methods is encouraged. Again it is emphasized that FMDF is not intended to be 
an alternative to DNS. The value in the FMDF approach lies in its applicability to more 
“practical” systems which are incapable to be treated by DNS. 


7.6 Validations via Laboratory Data 

In the experiments of Mungal and Dimotakis (1984), several runs involving different equiv- 
alence ratios were undertaken to elucidate the effect of this parameter on the product for- 
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Table 1: The computational times for 2D planar jet simulations. 


Simulation 

Grid resolution 

Normalized CPU timef 

Figure 

DNS 

1201 x 601 

242.5 

20, 21 

MlffHiiB 

201 x 101 

7.62 

19-21 

LES-FD 

201 x 101 

1 

20 


f Unit correspond to 760 seconds on a Cray-C90. 


Table 2: The computational times for 3D shear layer simulations. 


Simulation 

Grid resolution 

Normalized CPU timef 

Figure 

DNS 

217 x 289 x 133 

182.71 

22(b) 

FMDF 

55 x 73 x 34 

7.64 

22(b) 

LES-FD 

55 x 73 x 34 

1 

22(b) 


f Unit correspond to 655 seconds on a Cray-C90. 


mation. The authors define the equivalence ratio of the low speed freestream reactant con- 
centration Cq 2 to the high speed freestream value coj divided by the low speed to high speed 
stoichiometric ratio: 

c 02 /coi 


<t> = 


~T~ — C02/CQ1, 


(63) 


( C 02/0)1)3 

noting that the stoichiometric ratio is unity for the hydrogen-fluorine reaction. In all of the 
runs, the fluorine mole fraction was kept at 1% while the hydrogen mole fraction was varied 
at 1%, 2%, 4% and 8%. With the fluorine on the high speed side, this allows for equivalence 
ratios of 1, 2, 4 and 8. ‘Flip’ experiments were conducted in which runs with inverse values 
of the equivalence ratio were conducted (1, ^ and |). These were carried out by moving 

the 1% fluorine mixture to the low speed stream and the hydrogen mixture to the fast speed 
stream. 


Figure 25 displays contour plots of (a) the instantaneous temperature field and (b) the time 
averaged temperature field for a stoichiometric ratio of unity (1% fluorine on the high speed 
side and 1% hydrogen on the low speed side). In this simulation, the 8 species, 4 reaction 
finite rate chemical mechanism was utilized. Clearly evident in this figure is the presence 
of large scale coherent structures. Note that the peak temperature in the instantaneous 
field approaches, but is somewhat lower than the adiabatic flame temperature. Since the 
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chemistry is sufficiently fast, finite rate effects do not offer an explanation; rather the cause is 
explained by the nature of the filtering process. In contrast to the instantaneous temperature 
field, the peak values of the time averaged field are considerably lower than that of the 
adiabatic flame temperature, an intuitive fact noted previously by the authors as well as 
by others (Wallace, 1981). It is noted that a large number of individual particles (i.e. 
‘realizations’) do indeed approach the adiabatic flame temperature, a result reflected in Fig. 
26 which is a scatter plot of the paxticle temperature against the mixture fraction Z where 


z _ Slh + Cf 2°° ~ Cf 2 

c h 2 oo + c f 2oo 


(64) 


The presence of large scale structures in shear layers has been widely recognized (Brown 
and Roshko, 1974; Papamoschou and Roshko, 1986). Evidence of such structures in the 
present experiment is presented in the form of Schlieren photographs as well as temperature 
histories at several cross stream positions. Figure 27 is a corresponding time history of 
temperature at various cross stream locations calculated by the simulation. Each vertical 
increment represents temperatures ranging from ambient to the adiabatic flame temperature. 
Qualitatively, these time traces are strikingly similar to those generated in the experiment. 
One notable difference is that in the vicinity of the middle regions of the layer, there are points 
in time that the simulation exhibits near ambient temperature (cool fluid). While there is 
some evidence in the experiment, it is considerably more prevalent in the simulation. This 
is partially believed to be a consequence of the two dimensional nature of the simulation. 
Comparatively, the true three dimensional behavior exhibits smaller scales which tend to 
more effectively homogenize the fluid locally and prevent the presence of purely ambient fluid 
in this region. For this reason it is expected that the minimum values of the time averaged 
temperature in the vicinity y — 0 are slightly lower than those from the experiment. It is 
expected that a fully three dimensional FMDF simulation would alleviate this shortcoming. 
Another significant cause of this discrepancy noted by the authors is that the cold wire probes 
suffer from appreciable thermal lag and conduction error which tend to cause a ‘smoothing’ 
effect. 

In an effort to assess the capability of the simulation to quantitatively predict the experi- 
mental results, Fig. 28 shows plots of the time averaged temperature as a function of the 
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cross stream coordinate at the downstream measuring station. In this figure, the cross 
stream coordinate is normalized with the 1% temperature thickness, $ x , defined as the dis- 
tance between the points at which the cross stream mean temperature rise is 1% of the 
maximum mean temperature rise. The comparison between the experimentally measured 
values and those generated by the simulation procedure is favorable. Also in this figure is 
a plot of the mean temperature generated by the FMDF approach with the infinite rate 
reaction model. As expected, the infinite rate approach agrees quite well with the finite 
rate reaction mechanism. A considerable savings in computational time as well as memory 
is gained in the infinite rate approach, owing to the fact that only one passive scalar needs 
to be transported. Additionally, while 10 chemical time steps per hydrodynamic time step 
were required in the finite rate approach, the conserved scalar is not altered by the effects of 
reaction, eliminating the necessity for such expensive time splitting. As a result, while the 
finite rate code took 12, 700 cpu seconds on a Cray C-90, the infinite rate code used only 
1,990 cpu seconds on the same machine. With confidence in the infinite rate model, the 
remainder of the hydrogen-fluoride simulations were conducted using this approach. 

To demonstrate the ‘flip’ effect noted in the experiment, Fig. 29 shows the cross-stream 
temperature variation for (j> = 8 and <f) = |. Two primary observations are consistent with 
the findings of the experiment: (1) the peak value of the mean temperature differs although 
the adiabatic flame temperature is the same and (2) the peak value ‘shifts’ toward the lean 
reactant stream. Clearly, since the only difference between the two runs is the interchange 
of the low and high speed reactants, the mechanism for this observation lies in the fluid 
mechanics. An explanation is offered by Mungal and Dimotakis (1984) in the context of the 
Broadwell-Breidenthal theory (Broadwell and Breidenthal, 1982). This theory based upon 
the fact that the entrainment ratio differs from unity due to the non symmetric nature of 
the velocity profile (the entrainment ratio is defined as the ratio of the volume of high speed 
fluid to the volume of low speed fluid mixed within the layer). An additional explanation 
is presented based upon the PDF’s of Konrad (Konrad, 1977; Broadwell and Breidenthal, 
1982). In this approach, the flipping effect is attributed to the fact that the Konrad PDF’s 
have peaks near £/( 1 + £), where £ is the entrainment ratio. Regardless of the theoretical 
explanation, the important feature is that the FMDF simulation is capable of capturing the 
phenomena observed in the experiment. 
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Absolute temperature profiles for all equivalence ratios are presented in Fig. 30(a). These 
results exhibit characteristics similar to those presented in the experimental investigation. 
Most notable is the tendency of the flame to favor the lean reactant side as discussed above. 
Additionally, with the exception of the two cases for <j> = 1, the peak temperature is higher 
for equivalence ratios greater than one compared to the reciprocal equivalence ratio, a gen- 
eralization of the result discussed in the previous paragraph. Normalized mean temperature 
profiles are presented in Fig. 30(b). Consistent with the experiment, the peak normalized 
temperature reaches a maximum between the values of 1 and 2 for the equivalence ratio. 
These trends are more clearly portrayed in Fig. 31(a), which Mungal and Dimotakis (1984) 
refer to as ‘inferred’ temperature profiles. These plots are modified to reflect the temper- 
ature if the high speed reactant was fixed at 1 % molar concentration while the low speed 
stream was varied from through 8% to obtain the desired equivalence ratios. This figure 
clearly supports the conclusion drawn by the authors of the experiment that there exists an 
asymptotic limit to product formation as the high speed reactant is burned to completion. 
Similar behavior is exhibited in Fig. 31(b) in which inferred temperature profiles are shown 
for the situation in which the low stream reactant is fixed at 1% and the high speed reactant 
is varied to obtain the same equivalence ratios. 


Further assessment of the FMDF approach is made in Fig. 32(a) which shows the variation 
of the normalized product thickness with the equivalence ratio. In lieu of direct species 
measurements, product formation in the experimental investigation was inferred from the 
temperature measurements. While the simulation allows direct measurement of the species 
mass fractions, the product thickness as defined by Mungal and Dimotakis (1984) was used 
for consistency: 


& p \ = 
Sp2 = 


/ 

L 


+~ C r (T(y)) 

-oo cioAQ 

C r {T(v)) 

-oo C20AQ 


dy. 


dy. 


The distinction between these two definitions of product thickness lies in the free stream 
concentration used for normalization purposes (Mungal and Dimotakis, 1984). In the def- 
inition of S p i , the high speed mole fraction cq\ is used, while the low speed mole fraction 
C 02 is used in the normalization of S p2 . Figure 32(a) indicates the simulation captures the 
proper magnitude of product formation through all ranges of the equivalence ratio. At low 
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values of <f), the amount of product formed varies nearly linearly as the low speed reactant 
is consumed in reaction with excessive amounts of the high speed reactant. At high values 
of the equivalence ratio, the product thickness approaches an asymptotic value as reaction 
progress is inhibited by a lack of the high speed reactant relative to the amount of reactant 
in the low speed stream. This figure corresponds to the area under the curves in Fig. 31(a). 
Similarly, Fig. 32(b) shows the variation of the normalized product thickness with the 
inverse equivalence ratio. The corresponding temperature profiles are those of Fig. 31(b). 


8 Concluding Remarks 

The utility of the filtered mass density function (FMDF) and its application to LES of 
variable density chemically reacting flows is explored and demonstrated. The advantage 
of FMDF lies in its ability to resolve the correlations of the reactive scalars, density and 
temperature at the subgrid scale (SGS) level without resorting to closure models. The 
present formulation is developed for treatment of scalar variables only, thus modeling is 
required for correlations involving the velocity field (such as the Favre SGS stress and mass 
flux). This may be overcome by considering the joint velocity-scalar FMDF. In the present 
form the scalar FMDF methodology is capable of treating complex exothermic reacting flows 
and is powerful in a sense that it can be added to the existent CFD codes. 

The evolution equation for the FMDF is derived. Unclosed terms in this equation are identi- 
fied and models are used to represent their effects. The FMDF transport equation is solved 
numerically via a Lagrangian Monte Carlo scheme in which the solutions of the equiva- 
lent stochastic differential equations (SDEs) are obtained. These solutions preserve the Ito- 
Gikhman nature of the SDEs. The consistency of the FMDF approach, the convergence of its 
Monte Carlo solution and the performance of the closures employed in the FMDF transport 
equation are assessed by comparison with results obtained by direct numerical simulation 
(DNS) and by conventional LES procedures in which the first two SGS scalar moments are 
obtained by a finite difference method (LES-FD). These comparative assessments are con- 
ducted by implementations of all three schemes (FMDF, DNS and LES-FD) in temporally 
developing mixing layer and spatially developing planar jet under both non-reacting and 
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reacting conditions. In non-reacting flows, the Monte Carlo solution of the FMDF yields 
results similar to those via LES-FD. The advantage of the FMDF is demonstrated by its use 
in reacting flows. In the absence of a closure for the subgrid scalar fluctuations, the LES-FD 
results are significantly different from those based on filtered DNS. The FMDF predictions 
yield much closer agreement to the DNS. The FMDF methodology is also tested by com- 
parative assessment against the experimental data of Mungal and Dimotakis (1984) for a 
hydrogen-fluorine reacting mixing layer. It is shown that the FMDF method is able to accu- 
rately predict the temperature field. The “flip effect” and and other physical phenomenon 
associated with experimental results are also correctly predicted by FMDF. 

The results presented in this work demonstrate that FMDF provides a powerful tool for large 
eddy simulation of turbulent reacting flows. While the FMDF method is computationally 
more expensive than conventional LES method, it provides extensive flexibility for treating 
complex reacting flows over that of conventional LES methods without resorting to ad-hoc 
assumptions. Additionally, the FMDF approach can be applied to practical combustion 
systems where DNS is not feasible. Based on these observations, it is anticipated that LES 
of turbulent reacting flows with realistic chemical kinetics to be routinely conducted for 
practical applications in the near future. 
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Figure captions 


Figure 1. Cross-stream variation of the filtered density in temporally evolving mixing layer 
obtained by FMDF-1 at (a) t = 0, (b) t = 44. 

Figure 2. Vorticity thickness vs. time in the temporally evolving mixing layer. 

Figure 3. Contours of the filtered vorticity and temperature in the temporal mixing layer 
obtained by FMDF-1 (right hand side) and FMDF-2 (left hand side) at t = 44. The plots 
at the top represent the vorticity field and the bottom plots show the temperature field. 

Figure 4. Cross-stream variation of the mean filtered (a) pressure, (b) kinetic energy, (c) 
temperature, and (d) scalar at t — 44. 

Figure 5. Cross-stream variation of the percentage of the difference in pressure as obtained 
by FMDF-2 and FMDF-1 with different smoothing -in temporally evolving shear layer at 
t = 44. (I) Long-dashed line: no smoothing, NPG = 50 and A E = A; (I) Dotted-dashed 
line: smoothed with Gaussian filter, NPG = 50 and A E = A; (I) solid line: smoothed 
with box filter, NPG = 50 and A E = A; (I) dashed line: smoothed with local least square 
filter, NPG = 50 and A E — A; (I) Long-dashed thick line: smoothed with Gaussian filter, 
NPG = 200 and A# = A; (I) Dotted thick line: smoothed with Gaussian filter, NPG = 50 
and A e = 2A. 

Figure 6. Contours of the filtered values of the conserved scalar as obtained by (a) LES-FD, 
(b) FMDF at t = 44. 

Figure 7. Cross-stream variation of (a) filtered scalar, (b) generalized scalar variance in 
temporal shear layer at t = 44. 

Figure 8. Product thickness variation with time in the temporally evolving mixing layer. 

Figure 9. Cross-stream variation of the mean filtered density in 3D temporal shear layer at 
several spanwise locations. 

Figure 10. Scatter plots of filtered values of (a) the temperature, (b) the conserved scalar 
as obtained by M.C. and F.D. with FMDF-2 method in 3D temporal shear layer with 2D 
forcing. 

Figure 11. Cross-stream variation of the mean filtered (a) vorticity, and (b) pressure in 3D 
temporal shear layer at several spanwise locations. 

Figure 12. Scatter plots of filtered quantities as obtained by M.C. and F.D. in the 3D 
temporal shear layer with 3D forcing: (a) temperature obtained by FMDF-2, (b) conserved 
scalar obtained by FMDF-1. 

Figure 13. (a) Product thickness variation with time, (b) Cross-stream variation of the 
product mass fraction, in the temporally evolving mixing layer. 


47 



Figure 14. Cross-stream variation of the filtered temperature in temporally evolving mixing 
layer at t = 44. 

Figure 15. Vorticity thickness vs. time in the temporally evolving mixing layer. 

Figure 16. Contours of the filtered pressure in planar jet for different cases, (a) calculated 
by FMDF-1 with no smoothing of the filtered temperature, (b) calculated by FMDF-1 with 
smoothed filtered temperature with box filter, (c) calculated by FMDF-2 with no smoothing 
of the filtered temperature. 

Figure 17. Contours of the filtered temperature in planar jet for different cases, (a) calculated 
by FMDF-1 with no smoothing of the filtered temperature, (b) calculated by FMDF-1 with 
smoothed filtered temperature with box filter, (c) calculated by FMDF-2 with no smoothing 
of the filtered temperature. 

Figure 18. Contours of (a) particle number density and (b) fluid filtered density in planar 
jet for the case with uniform weights. 

Figure 19. Contours of (a) particle number density, (b) particle weighting and (c) fluid 
filtered density in planar jet for the case with variable weights. 

Figure 20. Temporal evolution of the product thickness, (b) Cross-stream variation of the 
product distribution at x/D = 14, in the spatially evolving reactive planar. 

Figure 21. Contours of the instantaneous subgrid scale unmixedness for the spatially evolving 
planar jet, (a) DNS and (b) FMDF. 

Figure 22. Product thickness variation with time in the 3D temporally evolving mixing layer, 
(a) low LES resolution, (b) high LES resolution. 

Figure 23. Cross-stream variation of the mean product mass fraction, (a) the mean values 
are obtained by averaging over streamwise and spanwise directions, (b) the mean values are 
obtained by averaging over streamwise direction and are plotted at two different spanwise 
locations. 

Figure 24. Temporal evolution of the vorticity thickness. The upper curves correspond to 
the cases with high LES resolution and the lower curves correspond to the cases with low 
LES resolution. 

Figure 25. Contour plots of (a) instantaneous Favre filtered temperature and (b) time 
averaged temperature fields for <j> = 1. 

Figure 26. Scatter plot of particle temperature versus mixture fraction for <f> = 1. 

Figure 27. Time history of instantaneous Favre filtered temperature at several cross stream 
locations. 

Figure 28. Cross stream variation of the normalized mean temperature for <j> = 1. 
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Figure 29. Cross stream variation of the normalized mean temperature for <j> = 8 and <f> = g. 

Figure 30. Cross stream variation of (a) absolute mean temperature and (b) normalized 
mean temperature for all equivalence ratios. 

Figure 31. Cross stream variation of the ‘inferred’ mean temperature profiles for (a) 1% high 
speed mole fraction and (b) 1% low speed mole fraction for all equivalence ratios. 

Figure 32. Normalized product thickness variation with equivalence ratio: (a) 8 P \ versus 
equivalence ratio and (b) 8 V 2 versus inverse equivalence ratio. 
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